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We present a novel, highly efficient yet accurate analytical approximation for the Green's function 
of a Holstein polaron. It is obtained by summing all the self-energy diagrams, but with each self- 
energy diagram averaged over the momenta of its free propagators. The result becomes exact for 
both zero bandwidth and for zero electron-phonon coupling, and is accurate everywhere in the 
parameter space. The resulting Green's function satisfies exactly the first six spectral weight sum 
rules. All higher sum rules are satisfied with great accuracy, becoming asymptotically exact for 
coupling both much larger and much smaller than the free particle bandwidth. Comparison with 
existing numerical data also confirms this accuracy. We use this approximation to analyze in detail 
the redistribution of the spectral weight as the coupling strength varies. 
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I. INTRODUCTION 

There is considerable interest to understand the ef- 
fects on the properties of a particle coming from interac- 
tions with an environment. Examples of such problems 
abound in condensed matter; the problem discussed here 
is that of an electron coupled to lattice vibrations, i.e. 
of electron-phonon coupling. For example, such coupling 
is believed to be relevant for understanding certain as- 
pects of the high-temperature superconductors' behavior 
in the underdoped limit (where the "particle" coupling 
to phonons is the doping hole already dressed by interac- 
tions with the electrons in the lower Hubbard band)ii2i^ 
but there are many examples of other materials char- 
acterized by strong-electron phonon coupling, including 
polymers like polyacetylene, nanotubes. Ceo molecules 
and other fuller enes^&i Other problems of the same gen- 
eral type regard electrons coupled to spin- waves of a mag- 
netically ordered system,^ or to orbitrons, for example in 
manganitesfSiiSi^ or to some combination thereof. 

In this work, we focus on the simplest Hamiltonian 
describing an electron on a lattice interacting with an 
optical phonon mode, namely the Holstein modelii^ 



k,q 

(1) 

The first term is the kinetic energy of the electron, with 
cj^ and Ck being the electron creation and annihilation 
operators. For the single dressed particle (known as po- 
laron, in this case) problem of interest to us, the spin of 
the particle is irrelevant and we suppress its index. 
is the free-particle dispersion. In all results shown here, 
we assume nearest-neighbor hopping on a d-dimensional 
simple cubic lattice of constant a and a total of N sites 
with periodic boundary conditions, so that 



£k = — 2t cos(/cia). 



(2) 



second term describes a branch of optical phonons of en- 
ergy (we set = 1 throughout this paper). 6^ and 
6q are the phonon creation and annihilation operators. 
The last term is the on-site linear electron-phonon cou- 
pling g c\ci{h\^hi)^ written in k-space. All sums over 
momenta are over the first Brillouin zone, — - < L < -. 

The quantity of interest to us is the Green's function 
of the single dressed particle, or polaron, defined as;i^ 



G(k,T) = 



'|I^[4(T)Ck(0)]|0), 



(3) 



where |0) is the ground state of the zero-particle system, 
which is the vacuum. T is the time ordering operator, 
and Ck(r) = e^'^^Cke"^'^^. Since H\^) = and Ck|0) = 0, 
Eq. simplifies to: 



G(k,r) = -z6(r)(0|cke- 
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where 0(r) is the Heaviside function. In other words, 
only the retarded part contributes. This Green's func- 
tion gives the amplitude of probability that an electron 
introduced in the system at r = and removed at a 
later time r, leaves the system in its ground state. The 
usefulness of this quantity can be appreciated using its 
Lehmann representatiomi^ 



G(k,u;)=^ 



|(a|4|0)|^ 
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but our results are valid for any other dispersion. The 



where and {E^} are the complete set of one- 

particle eigenstates and eigenenergies, 7Y|a) = Ea\a), 
Xlk^k^k|<^) = \c^). Thus, the poles of the Green's 
function give the whole one-particle spectrum^ while the 
associated residues, called quasi-particle {qp) weights, 
give partial information on the nature of the eigen- 
states. Moreover, the imaginary part of this Green's 
function, called the spectral weight, is directly mea- 
sured experimentally through angle-resolved photoemis- 
sion spectroscopy.^^ 

There has already been a large amount of work de- 
voted to understanding the properties of the Holstein 
polaron, and we briefly review some of it here. Most 
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of it is numerically-intensive work. Some examples are 
(i) exact diagonalization (ED) methods. ^^'^^ These are 
usually hampered by the fact that even for a finite 
lattice, the Hilbert space is infinite due to the infi- 
nite number of possible phonon configurations. Thus, 
some truncation scheme is needed, but for small phonon 
frequencies and/or large couplings, the CPU resources 
needed and run times become prohibitive. This has led 
to a number of (ii) proposals based on variational ap- 
proaches to decide which phonon configurations should 
be included, -^'^'-^^'-^^'^^'^-^'^^'^^'^^'^^ as weh as (iii) various 
Quantum Monte Carlo (QMC) methods. A brief review 
of these is provided in Ref. 27. Of special interest are the 
so-called diagrammatic Monte-Carlo simulationsS2iS2i2ft 
which calculate directly the Green's function in imagi- 
nary time, by numerically summing all diagrams in the 
perturbational expansion. We make extensive use of com- 
parisons with low-energy results of this method from 
Ref. ^3 While this method allows in principle the 
exact calculation of the Green's function, the require- 
ment of convergence for the propagator series usually 
means that only low-energy properties are shown. Fi- 
nally, there are methods suitable for some particular 
cases, such as density-matrix renormalization group for 
one-dimensional systems^^^ and dynamic mean-field the- 
ory (DMFT) for infinite-dimensional systems.^^ 

Of course, given the long history of this problem, many 
analytical techniques have been applied to it with varying 
degrees of success (for a review, see Ref. 33). First of all, 
the Green's function is known exactly in two asymptotic 
limits. If there is no coupling, ^ = 0, then the Green's 
function is that of the free electron. The ground-state is 
at 6o = —2dt and the spectrum consists of a continuous 
band extending from [—2dt^ 2dt] (for the tight-binding 
model). The so-called impurity limit, with t = 0, also has 
an exact solution, given by the Lang-Firsov formula^^ 

This can be viewed as the strong-coupling limit, since 
for t = 0, ^ becomes the important energy scale in the 
system. In this limit the electron is localized at one site 
in real space, therefore it is fully delocalized in k-space, 
and the Green's function is independent of k. The spec- 
trum has the ground-state (GS) at Eq = —g'^/^ and an 
infinite sequence of equidistant levels spaced by Q above 
it. This is extremely different from the free-particle spec- 
trum, and it is of considerable interest to understand not 
only how the ground-state evolves from —2dt to —g'^/Vt 
as the coupling g is increased, but also the evolution of 
all the higher-energy spectral weight from a continuous, 
finite- width band to an infinite set of discrete levels. Note 
that it is customary to describe the effective coupling as 
the ratio of the two asymptotic ground-state energies, 
using as a new parameter 



2dtVt' 



Most of the numerical methods reviewed above cal- 
culate only GS or low-energy properties, given the 
significant CPU time and numerical resources needed 
to calculate the whole spectrum. Very recently, 
several sets of whole-spectrum results have become 
available fS4iSii25i2fii21 however only for a few points in the 
parameter space, and generally for low dimensions. 

It is of obvious interest to find an analytical approxima- 
tion for the Green's function that is simple to estimate, so 
that the whole parameter space can be studied easily, but 
also with high accuracy. This is precisely what we pro- 
pose here (a short version of this work has been published 
in Ref. ,38. ). We call our approximation the momentum 
average (MA) approximation; its essence consists in an- 
alytically summing all the diagrams in the diagrammatic 
expansion, but with each diagram simplified in a certain 
way. Before introducing this method, we briefiy review 
here the other two simple (in terms of computational ef- 
fort) analytical approximations for the Green's function 
of the Holstein polaron, available in the literature. 

The first is the self-consistent Born approximation 
(SCBA), which consists of summing exactly only the non- 
crossed diagrams. The percentage of diagrams kept de- 
creases fast with increasing order (see Table P). If the 
coupling is small, the sum is dominated by the low order 
diagrams and SCBA works reasonably well. At strong 
coupling, the contribution of higher order diagrams be- 
comes essential, and SCBA is expected to fail (see below). 
In this approximation, the Green's function is written in 
terms of a self-energy: 

^ - ek - 2]scba(cc;) ^iri 
with the self-consistency condition 

SsCBA (^) = ]y X] GsCBA (k - q, - 11) . 

q 

Note that Sscba(^) is independent of k. This is a con- 
sequence of the simplicity of the Holstein model: if either 
the coupling g or the dispersion Vt were functions of the 
phonon momentum q, the SCBA self-energy would de- 
pend explicitly on k.^^ Sscba(^) can be expressed^^ as 
a function of the average of the free propagator over the 
Brillouin zone (BZ) and can be evaluated very efficiently. 

The other simple analytical approximation for the 
Green's function of the Holstein model is the general- 
ized Lang-Firsov (LF) expression.^^-^^ It is reminiscent 



order 


1 


2 


3 


4 


5 


6 


7 


8 


total 


1 


2 


10 


74 


706 


8162 


110410 


1708394 


SCBA 


1 


1 


2 


5 


14 


42 


132 


429 



TABLE I: Comparison between the total number of diagrams 
of a given order in the proper self-energy S(k, cj) vs. the 
number of diagrams of a given order kept within SCBA. 
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of the Lang-Firsov expression of Eq. 0: 
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This expression is exact in both asymptotic hmits A = 
{g = 0) and A ^ oo (t = 0), but less accurate for finite t 
and as we show below. 

The article is organized as follows. In Section II, we 
derive the momentum average approximation. Its dia- 
grammatic meaning is discussed in Section III, where 
we also estimate its corresponding spectral weight sum 
rules. We show that MA satisfies exactly the first 6 sum 
rules, but more importantly, it remains highly accurate 
for higher order sum rules. This is a strong argument 
in favor of its accuracy. The accuracy is gauged in more 
detail in Section IV, where we compare the MA predic- 
tions against the SCBA and generalized LF predictions, 
but also against a host of numerical results. This indeed 
demonstrates that the MA approximation is remarkably 
accurate for all parameter values, especially given its sim- 
plicity. In this section we also present some new results 
regarding various properties of the Holstein polaron. Fi- 
nally, Section V contains our summary and conclusions. 



II. CALCULATING THE GREEN'S FUNCTION 
A. Exact solution 

As is always the case for Green's functions, one can use 
the equation of motion technique to generate an infinite 
hierarchy of coupled equations for an infinite number of 
related Green's functions. We derive it here for the Hol- 
stein polaron. 

In the frequency domain, this approach is equivalent 
to using repeatedly Dyson's identity G{uj) = Go{uj) + 
G{u;)VGo{u;)^ which holds for any Green's operators 
G(cj) = [cj -n^ if]]-^, Gq(u) = [cj - Ho + ir]]-^ and 
for any Hamiltonians H = Hq ^V. As customary, we 
take V to be the electron-phonon interaction. Applying 
Dyson's identity once, we obtain: 



G(k,c.;) =Go(k,a;) 
where 



Go(k,w) = {0|ckGo(w)4|0) 



1 



(8) 



(9) 



is the free particle Green's function. We made use of the 
equality V'cj^lO) = Xlq <^k-q^ql^) defined a new 
Green's function: 

Fi(k,qi,c.) = (0|ckGHc[_q^6tjO). 

Fi is related to the amplitude of probability to start with 
the electron and a phonon at the initial time, and find 



only the electron in the system at the final time. Its own 
equation of motion relates back to G(k, uo) but also to 
a new Green's function with two phonons initially. In 
general, if we define: 

Fn(k, qi, . . . , q,, uo) = (0|ckG(c^)4_q^6^, ... |0), 

where qT = XI Qi is the total momentum of the n ini- 
tial phonons, using Dyson's identity we find its equation 
of motion to be (n > 1): 



Fn(k,qi, . . . ,qn,c^) = ^=Go(k - qT,cj - nQ) 
V N 



^Fn_i(k, qi, . . . ,qi_i,qi+i, . . . ,qn,Cc;) 



E 



+ ^n+i(k,qi,...,qn+i,tj) 



, (10) 



i.e. related to the Green's functions with n — 1 and n + 1 
initial phonons. Eqs. and (|Tnj) form the exact infi- 
nite hierarchy of coupled equations whose solution is the 
Holstein polaron Green's function G(k, cc;) = Fo(k, c^;). 

Obviously, this system of coupled equations can be 
solved trivially in the limit A = ^ = 0, in which case 
G(k, cj) = Go(k, cj) directly from Eq. 0. An exact so- 
lution equal to the Lang-Firsov result must also exist if 
t = 0. Indeed, in this limit all Green's functions become 
independent of all momenta, and Eqs. and sim- 
plify to: 



G{u) = Go{uj) l^gVNFi{uj 



Fn{u;) = gGo{u; - nQ) 



where Go{u;) = {co -\- ir])~^ . These recurrence equations 
can be solved in terms of continued fractions. We briefiy 
review the solution here. We suppress the functional 
notation and rewrite Fn = anFn-i + /^n^n+i, where 
an = ^Goioo - nQ), f3n = ^V^Go(cj - nQ). On phys- 
ical grounds we expect that F^+i becomes vanishingly 
small for a sufficiently large n, since it describes physical 
processes which are less and less likely. This allows one 
to solve these equations iterationally starting from this 
sufficiently large n: Fn ~ ani^n-i, to ffnd after solving 
for Fn-i, Fn-2, etc., that: 



-Fn 



(11) 



(^302 
1 - 



Allowing the continued fraction to be inffnite instead of 
truncated after n steps gives the exact solution. Recalling 
that Fo = G = Go(l + gy/NFi), one can now solve for 
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FIG. 1: The diagrammatic expansion for the self-energy. 



G. With the original notation, we find: 



1 



g'Go{uj)Go{uj-^) 
2g'^Go{uj -n)Go{uj -2n) 
1 



(12) 



1 



After some further work, this can indeed be shown to 
equal the Lang-Firsov expression of Eq. 0. 

This exact hierarchy of coupled equations [Eqs. 
and (|10|) ] can also be solved in the general case of fi- 
nite g and finite t by iteratively solving for Fi, then F2, 
etc., and removing them from this coupled system. It 
is straightforward to verify that the solution obtained in 
this case is the diagrammatic expansion, which can be 
partially resummed to the expected form: 



G(k,cj) 



1 



[Co(k, 



S(k,cj) 



(13) 



where the self-energy S(k, co) is the sum of all proper self- 
energy diagrams, the first few of which are shown in Fig. 
n While this solution as a sum of an infinite number 
of diagrams is exact, it is clearly not very useful if the 
sum cannot be performed. One typical strategy in such 
cases is to sum only a subset of these diagrams {e.g. the 
non-crossed ones, in the Self-Consistent Born Approxi- 
mation). This is reasonable when one can argue that the 
diagrams kept contribute much more than the neglected 
diagrams, which is not the case for SCBA in this prob- 
lem (see below). We propose a new strategy, explained 
in the next subsection, to find an approximative solution 
of these equations in the case of finite t and finite g. 



B. The Momentum Average Approximation 

To obtain an approximate solution in the case of finite 
t and g, we proceed as follows. We first note that G(k, cj) 
depends on Fi only through its average over the Brillouin 
zone /i (k, cj) = Fi (k, qi , cj) , 



G{k,uo)=Go{K^) [l+^V7V/i(k,c^)J . (14) 
We define the momentum averaged Green's functions: 

/n(k,cj) = ^ Fn(k,qi,...,qn,c^), (15) 

qi, ,qn 

where all momenta sums run over the BZ, and attempt 
to express Eqs. ([117)) in terms of these simpler quantities. 



by performing the corresponding momenta averages on 
both sides. While the first term on the right-hand side 
can be written in terms of fn-i exactly, the second term 
requires an approximation, which we choose to be: 

Fn+i(k,qi, . . . ,qn+i, 6^)6^0 (k - qT,c^ - nQ) 

^ Ngo{uj - nl])/n+i(k,cj), (16) 



qi,...,qri + i 



where 



^o(^) = ^^Go(k,cj) 



(17) 



is the free propagator momentum- aver aged over the BZ. 
One justification for replacing Go(k — qr,^ — nCl) by 
its momentum average go{^ — ^^) in Eq. ()16)) is that 
qT = qi takes, with equal probability, any value in 
the BZ. Moreover, in the impurity limit t = where all 
Green's functions are momentum independent, Eq. ([TB|) 
is exact. This suggests that our approach should be rea- 
sonable at least in the strong-coupling limit t <^ g. In 
a more practical sense, this approximation allows us to 
write /n(k, u;) in terms of /n-i(k, u;) and /n+i(k, u;) only, 
which is what we require to be able to obtain an analyt- 
ical expression for G(k, cj). We discuss the meaning and 
consequences of this approximation in more detail below. 
After the approximation of Eq. (|TH|) . Eqs. (|Tn|) be- 



come: 



/n(k,Cj) 



ggo{u; - nQ) 



N 



[n/n-i(k,cj) + A^/n+i(k,cj) 



Together with Eq. (|T^ these simplified recurrence rela- 
tions can be solved similarly to the t = case. We find: 



G(k,cj) = 



1 



(18) 



where the self-energy is, within the MA approximation, 

g'^Soi^ - ^) 



(^) 



1 



1 



(19) 



This is the main new result of this work. As discussed, 
the MA approximation becomes exact in the limit of zero 
hopping (t = 0), where ^o(^) Go{uj) = {uj -\- ir])~^ ^ but 
also for zero coupling, ^ = 0. In the following sections 
we show that the range of validity of this approximation 
extends well beyond these asymptotic limits, and that in 
fact the MA expression is reasonably accurate over the 
entire parameter space. 

Note that although the MA self-energy looks similar 
to the DMFT self-energy,^^ this is in fact a very different 
approximation. This issue is discussed in Appendix IXI 
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III. MEANING OF THE MA APPROXIMATION 

A. Diagrammatics 

To understand the diagrammatical meaning of the MA 
approximation we expand Eq. (|19|) in powers of g^: 

Sma(^) = 9^9o{uj - 1^) + / [2gl{u - n)go{co - 2n)] 
+ / [4gl{io - ^)gl{u - 2n) + 6gl{u - n) 

xgl{u - 2n)go{io - 3^)] + (20) 

and analyze the various terms. First, one can verify that 
the MA approximation generates the correct number of 
proper self-energy diagrams to all orders. Indeed, there 
is one term of order two terms of order 4+6=10 
terms of order and so on (see Table P). Moreover, 
to each of these terms we can associate an MA diagram. 
These have the same topology as the exact proper self- 
energy diagrams. The difference is that each free prop- 
agator Go(k, a;) in the exact self-energy diagrams is re- 
placed with a momentum averaged free propagator ^o(^) 
(with the correct frequency) in each MA diagram. 

Using Eq. (|T7)) . the first-order self-energy diagram is 
(see Fig. QJ: 



N 



^Go(k-q,cj-r2) = g'^go{uj -n), 



i.e. S^-*-^ (k, cj) = S^]^(a;), and thus MA is exact to first 
order. Differences appear from the second order dia- 
grams, where the two exact contributions (see Fig. QJ: 

Go(k-qi,cj-r2)Go(k-qi -^2,^-2^) 

qi,q2 

X [Go(k -qi.co-n)^ Go(k - q2,c^ - n)] (21) 
are replaced, within MA, by two equal contributions: 



2gY^{u;-n)go{^-2n)='^ ^ Go{quu;-Q) 

qi,q2,q3 

X Go(q2, - 2n)Go{qs,co - Q). (22) 



Comparing Eq. (|2j) to Eq. (|22|), we see that the MA self- 
energy diagrams have the correct number of free propaga- 
tors with the correct frequencies, however the momenta 
of the free propagators are un-correlated and individually 
averaged over. It is as if there is no connection between 
the momentum carried by a cloud phonon when it is emit- 
ted and when it is re-absorbed by the electron. Precisely 
the same holds for all higher order self-energy diagrams. 

To gain a better understanding of the difference be- 
tween the exact and the MA diagrams, let us further 
analyze the dependence on t of Eqs. and ((2^ . For 
t = the expressions are identical, because the free prop- 
agators become independent of momenta and, as already 



discussed, MA becomes exact. For finite t, we expand 
each free propagator as: 



Go(k,c^) = Go{uj) 1 + SkGo(c^) + (skGo(c^))' + 



where Go{u;) = {co -\- ir])~^ . Inserting this expansion into 
Eqs. (PT|) and and collecting powers of t yields the 
following. All 0{t) terms vanish in both the exact and 
the MA diagrams because they are proportional to a 
Xlq^q = 0- In fact, all odd-order powers in t vanish 
because Xlq^Tq^^^ = 0. Next, consider terms of order t^. 
Such terms arise either from expanding one free propaga- 
tor to 0{t^)^ or from expanding two different free propa- 
gators to 0{t). The former case leads to the same result 
for both the exact and the MA diagrams, and we obtain 
6 contributions proportional to = 2dt^ . The 

latter case, however, reveals a difference. Five of the six 
0(t^) such contributions from the exact diagrams van- 
ish because they involve propagators carrying different 
momenta, and, for example, Xlq^ q2 ^k-qi^^k-qi-qs = 0- 
The exception comes from the two outside free propaga- 
tors of the non-crossed diagram, which carry the same 
momentum and result into another ^k-qi ~ ^^^^ 

contribution. This is absent in the MA approximation, 
where different propagators always carry different mo- 
menta. It follows that the MA second order self-energy 
diagrams capture 6 out of the 7 finite 0{t^) contributions 
correctly. Similar considerations apply for higher order 
t powers and for higher order diagrams, differences be- 
tween the MA and the exact self-energy diagrams coming 
only from terms involving free propagators carrying equal 
momenta in non-crossed diagrams. However, the error 
from such missed terms becomes smaller and smaller as 
one goes to higher order diagrams because the percent- 
age of self-energy diagrams with one or more pairs of free 
propagators of equal momenta decreases exponentially. 

We conclude that the MA approximation captures 
most of the t dependence of each self-energy diagram, 
while summing over all diagrams. This analysis suggests 
that MA should be quite accurate for any finite g and t 
values. In the next section, we reinforce this conclusion 
by analyzing the sum rules of the spectral weight. 



B. Sum Rules 

For an even better idea of the accuracy of the MA 
approximation, we consider the sum rules for the spectral 
weight A(k,cj) = -;^ImG(k, cj): 



Mn(k) 



(23) 



For a problem of this type (single dressed particle), the 
sum rules can be calculated exactly to arbitrary order. 
The usual approach is based on the equation of motion 
technique. We review it briefly here in order to make 
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a few useful observations. The key step is to rewrite 



SO that we have: 



=0 



Mn(k) = -;^Im (i^ ) / dcje-^^^G(k,cj) 



r=0 



The integral is now simply G(k, r 0+). Using the 
definition of Eq. Q), we find, for any r > 0: 



7 \ n 

i— \ G(k, r) = -ie(r)(0|ckW"e-^«^c3,|0), 



so that the sum rules simphfy to: 

M„(k) = (0|ck:^"4|0). 



(24) 



These vacuum expectation values can be evaluated di- 
rectly with some effort. We find Mo(k) = 1, Mi(k) = ek, 
M2(k) = + g^^ M3(k) = 4 + 2g^e^, + ^2^, etc. 

One very important conclusion that can be drawn from 
this derivation is that these sum rules have the same func- 
tional dependence on the energy scales t^^l^g anywhere 
in the parameter space. Of course, in various asymptotic 
regimes, different terms dominate the overall value {e.g. 
M2(k) ^ g'^ if g:^t while M2(k) ^ el if g <^ t). How- 
ever, this shows that if one can evaluate the sum rules 
exactly in any asymptotic regime, for instance by using 
perturbation theory, the results hold true everywhere in 
the parameter space, even where perturbation fails. 

The second important conclusion one can draw from 
Eq. is that each term in Mn(k) is proportional to 

^PQnign-m-p ^ whcrc <p^m^n — TTi — p< n are integers, 
i.e. the sum rule Mn is a polynomial of total order n 
in the energy scales of the problem. More complicated 
dependence on t, g^ for example through exp(— ^^/l^^), 
simply cannot appear (see below). 

The first conclusion suggests an alternative derivation 
of the sum rules, which can also be used for the MA 
and SCBA sum rules. Namely, we use the diagrammatic 
perturbational expansion of the Green's function valid for 
g <^tto evaluate directly the integrals dco c<;^G(k, co) 
and retain the imaginary part. In this case, ^(k, co) = 
^n,i(k, cj), where I^n,i(k, uj), i = l,Sn are all 
Green's function diagrams of order n, i.e. containing n 
phonon lines. The multiplicity Sn = (2n — 1)!! = 1 • 
3 . . . (2n — 1). Each diagram /^^^^(k, cj) is a product of 
2n + 1 free propagators, summed over internal phonon 
momenta. Since for large frequency each Go(k, a;) 
{uj-\-ir])~^, it follows that for |cc;| oo, each I)n,i(k, uj) 
g'^'^/{uj -\- ^77)^^+^. Since any integrand that decreases 
faster than l/cj^ has a vanishing contribution to Eq. 
it follows that the diagrams of order n only contribute to 
the sum rules Mp(k) with p > 2n. Thus, even though 
G(k, Cfj) is the sum of an infinite number of diagrams, 
only a finite number of them, of low order, contribute to 
any given sum rule and the calculation can be done. The 
same holds true for the MA sum rules, the only difference 
being that the self-energy parts in the Green's functions 



diagrams are replaced with the corresponding MA self- 
energy parts. 

Let us analyze the differences between contributions of 
the exact and of the MA diagrams to the sum rules. It 
is straightforward to verify that: 



-— Im / 



2n+l 



n Go{q^^u;-n,)=g'^, 



-j poo 2n+l 

-Im / rfa;a;2n+i^2n -Q Go(q,,a>-fiO 

J — 00 

2n+l 



i=l 



and 



2n+l 



1 "'"^^ 
-Im / t/cjcj^^+V n ^o(qi,^-^i) 
^ -^-^ i=i 

>n+l 



i=l 



Both the exact and the MA diagrams of order n contain 
products of the general form 11?=^^ Go(qi, ^ — ^i). 
Some of the free propagators are actually Go(k, c^;) (al- 
ways the first and the last one, but there can also be 
intermediary ones connecting proper self-energy parts). 
All other free propagators have momenta dependent on 
the phonon momenta, which are summed over (in the ex- 
act diagrams), or are individually averaged over (in the 
MA diagrams). Since there is one-to-one correspondence 
between the number of exact vs. MA diagrams and their 
topologies, and since eq = 0, it follows that the exact 
and the MA diagrams of order n give precisely the same 
contributions to M2n(k) and M2n+i(k). Differences ap- 
pear in the contribution to M2n+2(k) if there is at least 
one pair of propagators in any of the self-energy parts 
of the exact diagram that carry the same momenta. In 
this case, the corresponding Cq.eq^. averages to 2dt^ when 
the sums over phonon momenta are carried for the exact 
diagrams, whereas these terms always average to zero 
for the MA diagrams. Since most free propagators in 
the self-energy parts have different momenta, such dif- 
ferences are quantitatively small. This is especially true 
for large phonon frequencies Vt^ where the contributions 
proportional to Vt captured correctly by MA scale like 
some power of 2n + 1. This analysis can be continued for 
higher sum rules, with similar conclusions. 

We can now summarize our findings. Only the 0^^ 
order Green's function diagram (the free propagator 
Go(k, cj)) contributes to Mo(k) and Mi(k). Since this 
is included correctly in the MA and the SCBA cases, 
both give the correct Mo(k) and Mi(k). In fact, this di- 
agram contributes an to Mn(k), which is always the 
leading power in t contribution. The 1^^ order diagram 
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FIG. 2: (color online) Ratio of MA (squares), respectively 
SCBA (circles) sum rules and the exact sum rules, vs. order 
n. Results are for ID and A; = 0, Q = 0.5t and A = 0.5, 1, 2. 
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FIG. 3: (color online) Ratio of MA (squares), respectively 
SCBA (circles) sum rules and the exact sum rules, vs. order 
n. Results are for ID and /c = tt, Q = 0.5t and A = 0.5, 1, 2. 



is also exact in both the MA and SCBA cases, there- 
fore both give the correct M2(k) and Ms{k.) sum rules 
as well. Differences appear from M4 onwards. Because 
SCBA only keeps 2 out the 3 second order Green's func- 
tion diagrams, and 5 out of 15 third order diagrams, etc., 
the leading terms in g are 2g'^ instead of 3^^ in M4(k), 
bg^ instead of 15g^ in M6(k), etc. This shows that SCBA 
fails all sum rules with n > 3 significantly in the strong 
coupling regime where the term proportional to g"^^ gives 
the most significant contribution to M2n(k) (similar con- 
clusions hold for odd sum rules). So even though SCBA 
always satisfies exactly the first 4 sum rules, it is a bad 
approximation for large g, where big discrepancies appear 
for n > 3. 

MA satisfies M^ik.) and M^(k.) exactly as well, be- 
cause these only depend on having the correct number 
and topology for the 2nd order diagrams. MA fails from 
M^ik.) onwards, however in a very different manner than 
SCBA. The leading term in g^ has the correct prefac- 
tor, because MA has the correct number of 3rd order 
diagrams. The error comes from the 2nd order diagram 
containing the non-crossed self-energy diagram, as dis- 
cussed. Indeed, instead of the exact sum rule: 

Me{k) =€i^g^ [54 + 6t^{2d^ - d) + AslQ + Sel^^ 

+ /(ISdt^ + 124 + 22£:k^ + 251]^) + 15/, (25) 

MA finds a sum rule M^^^ik) = Me{k) - 2dt^g^. The 
leading terms in the g ^ t and g ^ t limits are always 
exact (as expected, since MA becomes exact in these lim- 
its), and this is true for all orders n. For n > 6 some of the 
cross terms are missing, but these are a minority related 
to non-crossed diagrams, as explained. We therefore ex- 
pect the MA sum rules to remain highly accurate for 



larger n values. That this is indeed true for higher sum 
rules is shown numerically in Figs. [21and|31 where we plot 
the ratio of the MA respectively SCBA sum rules, and 
the exact sum rules of same order n. The results shown 
are for ID and /c = 0, tt, but similar trends are found in 
the other cases. For A: = all the spectral weight is at 
negative frequencies, therefore Mn{0) alternate signs for 
even/odd n, and this is reflected in the non-monotonic 
behavior with n. For /c = tt, most of the weight is at pos- 
itive frequencies and sum rules are always positive. The 
magnitude of the exact sum rules increases roughly ex- 
ponentially with n, for instance for A = 2 and ft = 0.5t, 
Mi4(0) = 119, 516, 000. For A: = and A = 0.5, both MA 
and SCBA are reasonably accurate, with a slight edge 
for MA at higher n. However, MA is clearly much more 
accurate for all the other cases shown, and its accuracy 
is expected to improve even more as one moves further 
into the asymptotic regions of weak or strong coupling. 

Before ending this section, one more issue needs to be 
addressed. It is obvious that the MA sum rules must cap- 
ture the contributions to Mn(k) proportional respectively 
to and g^ exactly, since MA is exact for both ^ = 
and t = 0. One may assume that this alone is sufficient 
for a good interpolation at finite t and g. That this is not 
so is shown by the generalized LF approximation, which 
is also exact for t = or ^ = 0. However, in this approxi- 
mation one finds M^^{k) = l,Mj^^(k) = Skexp (-^), 
etc. Ml and all higher sum rules have unacceptable de- 
pendence on the energy scales g and Q (see second obser- 
vation above), even though they become exact asymptot- 
ically. As shown in the next section, this approximation 
indeed performs rather poorly for finite t and g. 

We conclude that while MA satisfies exactly the first 
6 sum rules, it remains accurate for all higher order sum 
rules, and is asymptotically exact. This is another argu- 



8 



ment in favor of the accuracy of this approximation over 
the whole parameter space. In the next section, we com- 
pare the MA predictions to those of existing numerical 
simulations to further support this claim. 



IV. RESULTS 

We first list the explicit expressions of the momentum 
averaged Green's function ^o(^) of Eq. (|T7|) . For nearest- 
neighbor hoping, it is straightforward to derive: 



sgn(cj) 



(c) 



7r{uj ^ if]) \u -\- if] 



4t 



and 



1 



respectively, where 



4t 



CO -\- 2t cos{kza) + ir] ' 



and 



n7r/2 

Jo 



is the complete elliptical function of the first kindi^ 
These integrals can be performed numerically very ef- 
ficiently. More generally, for any free electron disper- 
sion ek to which corresponds the free electron density of 
states (DOS) po{e) = ^ Z)k^(^ ~ ^k), we have ^o(^) = 
J^oo depo{e){uj — e -\- ir])~^ (also see Appendix . 

The self-energy I]ma(^) is then calculated easily from 
Eq. by truncating the continued fraction to a high- 
enough level. For an error of order e it is necessary to 
go to a level with n such that ng'^go{uj — nQ)go{uj — (n + 
1)^) < e. Using the fact that for large enough n we can 
approximate ^o(^ — ^ ^o(^ — (^ + 1)^) ^ —1/(^^)7 
it follows that we must have 



n > 



I 9' 



This result is expected, since ^ is roughly the average 
number of phonons in the polaron cloud (see below). This 
condition shows that all diagrams with at least that many 
phonons have to be included. In practice, we always use n 
large enough so that the change in I;ma(^) after doubling 
n is below a threshold much smaller than rj. All MA error 
bars in the figures we show are less than the thickness of 
lines/symbols used for the plots. 

With an explicit form for Ema(^) we are now in a 
position to calculate the Green's function GmaO^^ ^) and 
extract various polaron properties. 



A. Polaron Ground State Properties 

We begin by discussing polaron ground-state proper- 
ties. Most of these are already known from numerical 
studies, but they give us an opportunity to further test 
the accuracy of the MA approximation. Given the sim- 
plicity and efficiency of the MA approximation, we can 
also investigate higher dimensionality and larger parame- 
ter ranges than typical numerically intensive approaches. 
In this section, we use for comparison ID and 2D numeri- 
cal results obtained from diagrammatic Quantum Monte 
Carlo (QMC) simulations.^^ 

For /c = 0, we track the energy and weight of the lowest 
pole in the spectral weight, which give the ground- state 
energy ^0 and the ground-state quasi-particle weight 
Zq = |(G5'|cj^^o|0)p. Using the Hellmann-Feynman 
theorem, we then find the average number of phonons 
in the ground-state to be 



Nph = {GS\Y,blb^\GS) 



dEo 

on ■ 



(26) 



Note that one can also calculate the correlation function: 



{GSlY^clc, (bl + bi)\GS) 



dEo 
dg 



(27) 



just as easily. We do not show it here because we do 
not have the corresponding numerical data for the com- 
parison, however some typical results for this quantity 
are shown at the end of this section. Also note that 
all these quantities can be calculated similarly for other 
eigenstates. We will show such results in other sections. 

We also show the effective mass, m*. Because the MA 
self-energy is momentum independent for this simple Hol- 
stein model, one has>2£ 



m 
m 



1 



du 



(28) 



This result also gives us a consistency check on our calcu- 
lations. For MA we generally show effective mass results 
obtained from the first equality. 

A comparison of the results for these four quantities 
as obtained with QMC (black circles) and the differ- 
ent approximations is shown in Fig. ^ As expected, 
SCBA (blue line) fares well at small couplings A but very 
poorly at strong couplings. The generalized LF (green 
line) is asymptotically exact, but quite wrong at finite 
A. Because these are GS properties, we can also use 
perturbation theory in the two asymptotic limits to es- 
timate them. At weak coupling we use the Rayleigh - 
Schrodinger perturbation theory (RS, violet line) which 
gives the lowest energy for a state with total momentum 
kasi^ 
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(29) 
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FIG. 4: (color online) Ground state results in ID. Shown as 
a function of the coupling A are the ground state energy Eq; 
the qp weight Zo; the average number of phonons A^ph and 
the effective mass m* on a logarithmic scale. The left panels 
correspond to 0/t = 0.1 and the right ones to 0/t = 0.5. 



FIG. 5: (color online) Ground state results in 2D. Shown as 
a function of the coupling A are the ground state energy Eq; 
the qp weight Zq; the average number of phonons A^ph and 
the effective mass m* on a logarithmic scale. The left panels 
correspond to 0/t = 0.1 and the right ones to Q/t — 0.5. 



At strong couplings we use the second order perturbation 
theory result (PT, cyan line):iSi^ 



- eke ^ 



(30) 



^ 9^ 
The OS energy is simply Eq = E\^=o. A/'ph is obtained 



from £^0 as before; m* can be evaluated from the second 
derivative of E\^ with respect to k, and Zq is extracted 
from the effective mass [Eq. (pS)) ]. 

Fig. shows that one or the other of these perturba- 
tional values describe the OS energy Eq quite well, es- 
pecially for the larger Q value. However, the agreement 
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FIG. 6: (color online) Ground state energy ^o; the qp weight pjQ 7. (^olor onhne) Ground state energy ^0; the qp weight 
Zo ; the average number of phonons N^i, and the effective mass ; the average number of phonons A^ph and the effective mass 

m* as a function of A and of Q/t, for d = 1. m* as a function of A and of Q/t, for d = 2. 



for the other quantities is somewhat poorer, especially 
at strong couplings (PT and LF give identical results for 
Zq). Part of the reason is that the term in Eq. ((5n|) has 
in fact a more complicated dependence of g and 17, which 
is only asymptotically equal to the one used here^ More 
importantly, neither perturbational theory describes well 
the transition area, or can be easily applied to high en- 
ergy states. 



Clearly, MA (red line with red symbols, for easier com- 
parison with QMC data) has the best agreement with 
the QMC data. As expected from the sum rule analy- 
sis and the discussion on the convergence of Uma, MA 
improves for larger 17. The worst disagreements we ever 
found are the ones shown in the intermediary A regime 
for Vt/t = 0.1. Even there, the error in the GS energy 
is always below 5%. The qp weight has a more signifi- 
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cant disagreement, however note that it indeed becomes 
asymptotically correct for A < 0.2 and A > 0.8. The 
second claim is supported by the m* data, which in- 
deed shows convergence towards the expected PT values. 
Most importantly, even though it is quantitatively some- 
what wrong in this intermediary regime for small the 
MA approximation is the only one that clearly captures 
the crossover from the large to the small polaron, which 
is accompanied by the collapse of the qp weight and the 
increase in the number of trapped phonons and thus of 
the effective mass. 

The agreement with the QMC data is significantly im- 
proved for 2D polarons, as shown in Fig. Here MA 
gives excellent agreement at all couplings A. This is all 
the more remarkable when considering what a numeri- 
cally trivial task it is to evaluate the MA results com- 
pared to the QMC simulations. The physics is similar to 
that seen in ID, however the cross-over from the large 
(light) polaron at weak couplings, to the small (heavy) 
polaron at strong couplings, becomes somewhat sharper, 
especially for smaller Q/t values. 

Given the simplicity of MA, we can also generate con- 
tour plots of these quantities as a function of both A and 
Q/t, and investigate the entire parameter space. Such 
results are shown in Figs. El and d for d = 1 and 2, re- 
spectively. The cross-over from the large to the small po- 
laron can now be tracked (for instance, from the collapse 
of Zq) and one can quantitatively trust the results to a 
high degree. The transition occurs for A ~ 1, although 
this shifts to lower values as Q decreases. 

We also show MA results in 3D, see Fig. |H1 Here we do 
not have QMC results for comparison, however the good 
agreement with one or the other perturbational theories 
for most coupling strengths suggests that the MA results 
are probably even more accurate than in 2D. This is con- 
sistent with the expectation of improved agreement in 
higher dimensions. The crossover from large to small po- 
laron becomes even sharper, especially for lower Cl. It is 
still located in the neighborhood of A 1. 

To the best of our knowledge, the highly accurate 
three-dimensional results shown in Fig. |S1 are the first 
of their kind. This is likely due to the fact that the nu- 
merically intensive techniques require far too much com- 
putational effort to investigate such cases ^2ii2i22i^ 

Using the Hellmann-Feynman theorem like in Eqs. 
((^ and (P7)) also allows us to separate the individual 
contributions of 
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FIG. 8: (color online) Ground state results in 3D. Shown as 
a function of the coupling A are the ground state energy ^o; 
the qp weight Zq; the average number of phonons A^ph and 
the effective mass m* on a logarithmic scale. The left panels 
correspond to Q/t = 0.1 and the right ones to Q/t = 0.5. 



and 



{GS\g^ci 



h]\GS) 



to the total OS energy. Plots of these individual contri- 
butions as a function of coupling strength A are shown 
in Fig. ini for d = 1,2,3. As expected, the kinetic en- 
ergy is close to —2dt at weak couplings, but it becomes 
vanishingly small in the strong coupling limit, where the 
polaron becomes very heavy. The phonon energy E^\^ 
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FIG. 9: (color online) GS expectation values for the electron 
kinetic energy (violet), the phonon energy (cyan), and the 
electron-phonon interaction (orange), as a function of A, for 
Q 0.5t and (i= 1,2, 3. 



increases roughly like g"^/^ in the strong coupling limit, 
where A^ph ^ g'^ /Vfi. It follows that the decrease in the 
total GS energy is due to the interaction term, as ex- 
pected. Note that this energy is proportional to the cor- 
relation of Eq. ((TT)) . Since Eq ^ —g'^/^ in the strong 
coupling limit (see agreement with PT results), it fol- 
lows that this correlation becomes asymptotically equal 
to —2g/Cl in the strong coupling limit. 



B. Low energy states: momentum dependence 

We can also calculate the same properties for the lowest 
energy state corresponding to each given momentum k 7^ 
0, to find the low-energy behavior of the polarons. In 
this section we present comparisons with available QMC 
resultsi^ for ID and 2D systems. 
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FIG. 10: (color online) Lowest polaron eigenenergy for a 
given /c, E^^ and the corresponding qp weight Zk and av- 
erage phonon number N^\^{k). Results are for ID, VLji — 0.5, 
A = 0.25 and 1 (for E]<^ and Z^) as well as 1.96, for A^ph(/c). 
Only half of the BZ is shown. 



In Fig. ^Jwe show ID results for the polaron disper- 
sion £^/e, the associated qp weight Z/^ and average phonon 
number for two couplings. For the very weak 

A = 0.25, we see that MA and SCBA are equally good at 
small /c, however for large k MA overestimates the energy, 
such that the continuum which is expected to appear at a 
distance Vt above the GS energy £^0, is in this case pushed 
somewhat higher. We will return to a discussion of this 
discrepancy later on. As expected, the qp weight is large 
for small /c, where the main contribution to the eigen- 
state comes from the free electron state c^|0). However, 
Z/e goes to zero for larger /c, since these are primarily lin- 
ear combinations of states of type 6j|0), as confirmed 
also by the the average phonon number of about unity. 
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FIG. 11: (color online) Lowest polaron energy for a given k, 
^k, and the corresponding qp weight Zk and average phonon 
number A^ph(k). Results are for 2D, Q/t = 0.5, A = 0.845. 
Three high- symmetry cuts in the Brillouin zone are shown. 



Note that the RS perturbation also works well for small 
/c, as expected. It however breaks done at a finite k where 
Ck ~ —2t-\-Q, i.e. the free electron dispersion crosses into 
the continuum of electron plus one phonon states. Here, 
RS predicts an unphysical peak in the polaron dispersion 
(see the denominator of ((221)) and it fails for larger k. 

The second coupling A = 1 is roughly in the cross-over 
regime, see also Fig. ^ Here MA gives a much better 
agreement with QMC than SCBA or the perturbational 
theories. The polaron bandwidth is already renormalized 
and well below the weak coupling value of In fact, 
as we show later, there is another bound state between 
these states and the continuum. This is not captured in 
SCBA, which always predicts a polaron bandwidth of Q 
with a continuum above, and roughly between zero and 
one average number of phonons as k increases from to 
TT. For the strong coupling A = 1.96, low-energy proper- 
ties become almost /c-independent, as expected since the 
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FIG. 12: (color online) (a) Energy gap A = — between 
GS and first excited /c = state. A second bound state is 
stable if A < Q (here Q = 0.5t); (b) Line below which a second 
bound state appears: ED^^ (circles) and MA (squares); (c) 
MA qp weight of the second bound state when stable (red 
squares), and that of the GS (black squares), for Q/t = 0.5. 
The circle is the one QMC result available for the qp weight 
of the second bound state. These results are for ID. 



Lang-Firsov impurity limit is being approached. 

A second such comparison is possible for the 2D case 
with Q = 0.5t and A = 0.845, where QMC data is avail- 
able, see Fig. ^2 This coupling is on the weak side of 
the crossover, with the GS qp weight still large, Zq ^ 0.6. 
Here MA is already doing better than SCBA. As in the 
ID weak-coupling case, one can again see that the MA 
polaron bandwidth is slightly larger than Q. MA also 
somewhat overestimates the average number of phonons 
for large k values. Overall, given that all this data is 
in the crossover regime where the MA is at its worst, 
one can conclude that MA is also reasonably accurate in 
capturing low-energy polaron behavior. 



C. High-energy states 

The main motivation in trying to find an approxima- 
tion for the Green's function is that this quantity gives 
not only low-energy information, but the whole spec- 
trum. Here we compare MA predictions with various 
high-energy results available in the literature. Unfortu- 
nately there are much fewer of these, because the com- 
putational effort to obtain the whole spectrum through 
the usual numerical approaches is generally forbidding. 

We begin with a comparison against exact diagonaliza- 
tion (ED) ID data, from Ref. [l3|. The results are shown 
in Fig. El As already discussed, for weak coupling there 
is a continuum starting at + but for stronger cou- 
pling a second so-called bound state appears below the 
continuum. In Fig. IT^ a) we track the energy Ei of the 
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second k = state. For small couplings, the data ac- 
tually shows the maximum DOS in the continuum, not 
its edge (the maximum is generally located close to the 
lower edge. This data shows again that MA somewhat 
overestimates this energy, which should be ^ Q). When 
El < Eq -\- 1], there is a true discrete state. Note that 
panel (a) is in very good quantitative agreement with 
similar data shown in Fig. 8 of Ref. 17. The only dif- 
ference is for strong coupling, where the ED data shows 
El > Eq -\- Q again, however, with significant finite-size 
dependence on the chosen Hilbert space cutoff. 

We can thus find the coupling g/t where the second 
bound state appears, for different values of Q/t. This 
line is shown in panel (b), together with the ED results. 
The agreement between the two data sets is excellent, 
even at small Q/t values where we expect MA to be less 
accurate. We also show in panel (c) the qp weight of this 
second bound state, where stable. This data is not given 
in Ref. IT^ however one QMC point is available in Ref. 
I30I in good agreement with the MA prediction. 

We now move to comparisons for the entire spectral 
weight A{k^u;). In Figs. IT!^ IT^ and IT^ we show compar- 
isons for a ID system with Q = OAt and three different 
coupling strengths A = 0.5, 1 and 2, respectively. In each 
case, data for 5 values of /c, namely 0, 5, and tt 

are shown. The numerical data (black line) is obtained 
using a variational method by De Filippis et aL^^ Nu- 
merical data obtained through exact diagonalization of 
a finite system and by QMC, for the same parameters 
but somewhat different k values, is also presented by Ho- 
henadler et al. in Refs. IsslIsS These sets of numerical 
data are in good agreement with one another. 

In all three cases the agreement between MA results 
and the numerical data is very good. As expected, it is 
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FIG. 14: (color online) ID spectral weight A{k,(jj) vs. cj, 
for k — 0,^f,-^ and tt. MA results (red line) vs. data 
from Ref. |^ (black line). Parameters are Q = 0.4t, A = 1, 

77 = o.ia 



best for the largest A, but even for the smaller A values, 
which are just below and within the cross-over region, the 
agreement is very satisfactory. For A = 0.5 and /c = 
(upper panel of Fig. ITSj) we see the polaron state as a 
Lorentzian peak (a broadening rj = {)AVt was used) which 
accounts for most of the weight, and a small continuum 
at a higher energy. MA overestimates the gap between 
the two, which should be Vt. As k increases, the polaron 
peak disperses but also looses significant weight, as dis- 
cussed in the previous section. Most of the weight is now 
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FIG. 13: (color online) ID spectral weight A{k,oo) vs. 00, for 
k — 0, f , ^- results (red line) vs. data from 

Ref. O (black line). Parameters are Q = 0.4t, A = 0.5, 
T] = O.IQ. 



FIG. 15: (color online) ID spectral weight A(k^ijj) vs. u;, 
for k = 0, ^ , f , ^ and tt. MA results (red line) vs. data 
from Ref. |2J (black line). Parameters are — 0.4t, A = 2, 

T] = o.ia 
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FIG. 16: (color online) MA ID spectral weight A{k = 0,uj) 
vs. 00, for Q = 0.4t, A = 1, and r]/Q = 0.1,0.05,0.02,0.01. 
The first three peaks are discrete states (lorentzians) whereas 
the fourth marks the band-edge singularity of the continuum. 



in the high energy continuum, located roughly near the 
corresponding value. This simply shows that these 
higher energy states are not significantly affected by this 
rather weak coupling. The VM data shows somewhat 
more structure in these continua than the MA data, but 
most of the weight occupies similar frequency ranges. 

For A = 1 and k = (upper panel of Fig. EI), the MA 
data shows 3 Lorentzian peaks plus a continuum start- 
ing at u/t = —1.6. For the rather large 77 used it is hard 
to distinguish which peaks come from individual poles, 
and which are true continua. This can be easily done by 
studying their behavior as a function of the broadening 77, 
as shown in Fig. ^1 The height of peaks corresponding 
to discrete states scales precisely like I/77, as expected 
for lorenzians. The continuum is affected very little by 
changes in 77, except the peak near its lower edge where 
the finite r] smoothes out a singularity in the DOS. Since 
this singularity is not of the l/co type, its scaling with r] 
is different from that of the Lorentzians. The two lower 
states are closer to one another than f^, however the MA 
data shows no sign of the continuum that is expected 
to start at £^0 + ^- Note that the numerical data in Fig. 
ICT a) shows more structure, that could be consistent with 
this continuum. We will address the issue of this contin- 
uum below. As k is increased (see Fig. [T^ the low-energy 
peaks show some dispersion, but with a strongly renor- 
malized bandwidth. At higher k most weight shifts again 
at high energies, in a rather broad continuum. Finally, 
for A = 2 and /c = 0, Fig. [T5I shows even more discrete 
peaks spaced by ft. The GS is at £^0 ^ — 4.25t, but its 
weight is so small that it cannot be seen on this scale, 
unless T] is decreased significantly. A continuum is seen 
above co = — 1.6t. As k is now increased, there is almost 
no dispersion of the discrete peaks, however the weight 
shifts again to an even broader high-energy continuum. 

The issue of the continuum at Eo-\-Q in the exact case, 
and of its absence in the MA approximation for moder- 
ate and large couplings can be understood from Fig. El 




FIG. 17: (color online) 2D spectral weight A(k = 0,cj) 
for Q = 0.5t, T] = O.Olt and various A values, from exact 
diagonalization^'^ (black) and MA (red). For A > 1.125, the 
arrows point the GS location. The insets show the same data 
on a smaller scale, so that low- weight states are more visible. 
The arrows show the GS and the state appearing at precisely 
Q above GS, in the ED results. 



Here we show a comparison of A{k = 0, 00) in 2D ob- 
tained from exact diagonalization^ vs. MA results, for 
various couplings A. The first remark is that MA cap- 
tures quite well all the large- weight features, both as far 
as their energy and their weight are concerned. This is 
expected, given the good sum rules agreement demon- 
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strated previously. However, the ED results clearly show 
more states than MA predicts. There is always a low- 
weight peak at precisely Eq -\- Q (the GS and this state 
are marked by arrows in the insets). For couplings of 
up to A ~ 1, this peak is followed by several nearby 
peaks with comparably low- weight, which can be argued 
to be part of the expected continuum. For larger A, how- 
ever, only the state at + ^ can still be seen, although 
more states suggesting more continua are seen between 
the large higher-energy peaks. The gradual disappear- 
ance of the first continuum is not surprising, since one 
expects its width to narrow exponentially as the cou- 
pling increases. Moreover, one expects that the largest 
contribution to this continuum is from states with one or 
more phonons, explaining their low qp weight. 

As far as MA is concerned. Fig. [T7I suggests that for 
couplings where there is more than one discrete state, the 
very little weight in the + ^ and similar higher-energy 
continua is absorbed in the discrete states predicted by 
MA. This is consistent with the systematic up-shift of 
the MA peaks compared to the ED data. 

In fact, it is straightfoward to see that the MA ap- 
proximation can only predict a continuum starting at 
—2dt -j-Q. A continuum is signaled by a finite imaginary 
part of Hma('^), and the lowest frequency where this can 
occur is that for which go{u; — fl) acquires a finite imag- 
inary part [see Eq. CHI)]- However, the imaginary part 
of ^o(^) is proportional to the total density of states of 
the free model, i.e. it is finite for co G [—2dt^ 2dt] for 
nearest-neighbor hopping. It follows that the MA con- 
tinuum always starts precisely at —2dt-\-Q. This explains 
why for small A, where there is only one peak below this 
continuum, the gap between the two is somewhat larger 
that the expected Q value: the GS energy decreases be- 
low —2dt with increasing A, whereas the continuum edge 
is pinned at —2dt + 1], in the MA approximation. As the 
coupling increases, bound states start to split from this 
continuum, and spectral weight is shifted to lower en- 
ergies, in good agreement with the sum-rule predictions 
of the exact solution. These new bound states have to 
account for the (small) weight that is present in lower en- 
ergy continua, in the exact solution, and this is precisely 
what Fig. El shows. Clearly, a self-energy that would 
account for these continua as well would have to be a lot 
more complicated than that of Eq. (|19|) . 

This shows again that MA is remarkably successful 
in predicting the main features of the Green's function, 
given its simplicity and trivial numerical cost. This 
should make it (and generalizations of it to other models) 
of large interest for comparison against experiments. 

In the following, we use MA to investigate more prop- 
erties of the Green's function. To our knowledge, there 
are few or no similar results with similar accuracy pub- 
lished in the literature. We begin with Fig. where we 
plot the qp weight and average phonon numbers for a few 
of the higher-energy peaks, once they appear below the 
continuum. For comparison, we also show the already 
discussed GS results (black line). Results in higher di- 




FIG. 18: (color online) qp weight and average number of 
phonons in the GS (black line) and the next three higher 
/c = bound states, when they become stable according to 
MA. Results are for ID, Q = 0.5t. 



mension are qualitatively similar and we do not show 
them here. Unlike for the GS, both these quantities are 
non-monotonic functions of A for all higher-energy bound 
states. Each of these states disperses with /c, like in Figs. 
[T^ and [TT)1 f more data for this is shown below), so an ef- 
fective mass can be associated with each such band. This 
effective mass satisfies m*/m = 1/Z, and therefore also 
shows non-monotonic behavior, first decreasing and then 
increasing as A is increased. The average phonon number 
in the nth state must approach g'^ /Vfi^n asymptotically, 
as can be verified in the Lang-Firsov limit. This is in- 
deed observed in Fig. however the plateaus seen at 
moderate A suggest some cross-over from one to another 
type of wavefunction associated with these higher levels. 
As A ^ oo, an infinite sequence of such bound states 
appear, as expected in the Lang-Firsov limit. 

For a better illustration of the appearance of these 
bound states and of their evolution, we show contour 
plots of the spectral weight A(k, cj). We begin by plot- 
ting A(k, uo) as a function of k and cj, for fixed parameters 
g^t and Vt. In Fig. ^Jwe show ID results correspond- 
ing io Vt = t and A = 0.4, 1 and 2, respectively. Only 
half of the BZ is shown, since time-invariance guaran- 
tees that G(k, cc;) = G(— k, cj). Each of these MA con- 
tour plots takes below ten seconds to generate. Note 
that similar plots for the same parameters were provided 
by Hohenadler et al. in Ref. |3^ based on a cluster 
perturbation theory approach. The agreement between 
the main features of our and their plots is excellent. As 
expected, their data does show a few more low-weight 
features at lower energies, below —2t + (7 = — t, in this 
case, where our continuum starts. Such contour plots 
are richer versions of plots like those shown in Figs. 
El and El They illustrate basically the same points, 
although one can now also see clearly the dispersion of 
various features. For the lowest A, the free electron dis- 
persion Ck = —2tcos{ka) is still almost visible, except 
that electron-phonon interactions split it into the lower 
polaron band, and the higher continuum. This contin- 
uum is not featureless, instead one can already see weight 
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FIG. 19: (color online) Contour plots of the ID spectral 
weight A{k,co) as a function of k and co. The intensity scales 
are shown to the right of each plot. Parameters are Q — t, 
A = 0.4, 1 and 2, respectively, and 77 = 0.02. These results are 
in excellent agreement with the results of Hohenadler et al, 
Ref . and the GS results of Bonca et al. , Ref . il7i . 



accumulating near its lower edge. As A increases, a new 
bound state will split off from this. Indeed, this is seen 
for A = 1, where there are 2 bound state below the con- 
tinuum starting at —1 (for these parameters). The band- 
width of each of these states is now narrowed below Vt. 
The weight in the continuum at higher energies is re- 
distributed suggesting the impending formation of yet 
more bound states. Indeed, the A = 2 data shows 4, 
even narrower bound states below the continuum, which 
is showing signs of further fragmentation at multiples of 




FIG. 20: (color online) Contour plots of the 2D spectral 
weight A{k^ijj) as a function of k and a;, along several cuts 
in the BZ. The intensity scales are shown to the right of each 
plot. Parameters are Q = 2t, A = 0.5, 0.945 and 2, respec- 
tively, and Tf — 0.02. The middle panel is in excellent agree- 
ment with the results of Hohenadler et al. given in Ref. l35l . 



the phonon frequency. 

Similar behavior is expected, and indeed seen, in higher 
dimensions. Here we only show similar 2D contour plots, 
for 1] = 2t and A = 0.5,0.945 and 2, in Fig. \M The 
middle pan el again agrees very well with data shown 
in Ref. Issl In this case, the MA continuum starts at 
— 4t + 1] = —2. As A increases, we see again first 1, then 
2 and then 4 bound states below the continuum. Their 
bandwidths decrease with increasing A, so that the low- 
est band for A = 2 is already almost dispersionless, even 
though its weight still varies with k. As in the ID case, as 
the interaction becomes stronger, the weight in the con- 
tinuum also redistributes itself, with strong resonances 
seen around multiples of the phonon frequency. 
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since within the MA approximation this is given by: 
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FIG. 21: (color online) ID results for = 0.5t and ry = 0.04t. 

(a) A{k = 0,0;) vs. uj and A; (b) A{k = 7r,u;) vs. a; and 
A, on a linear scale; (c) A{k — 7r,a;) vs. and A, on a 
logarithmic scale; (d) total spectral weight A{ix}) vs. a; and A. 
The intensity scales are shown to the right of each plot. 



Another way to understand the dependence on the cou- 
pling A (or any other parameter) is to plot a contour of 
the spectral weight A(k, uo) vs. A and uo for a fixed value 
of k. Such a task is equally trivial at the MA level. In 
fact, one can also just as easily calculate and plot the 
total density of states, or spectral weight: 



A{w) 



In Fig. we show 4 such contour plots for the ID 
case. The uppermost one shows A{k = ^^uo) vs. uo. For 
A = (noninteracting case), only one state exists at — 2t, 
as expected. As the coupling turns on, the energy of 
this state (the ground-state) decreases, but /c = weight 
is also transferred to higher energies, due to hybridiza- 
tion with the states in the electron-plus-one-phonon con- 
tinuum. The MA continuum here starts at — 1.5t. For 
moderate and larger A one can clearly see how weight 
is re-arranged inside the continuum as A increases, and 
new bound states split from it and move towards lower 
energies. The apparent "break" in the slope of the GS 
energy, as A increases, is now seen to occur when the first 
bound state approaches the GS, and is suggestive of an 
avoided crossing. From this point on the GS lowers its 
energy much faster, but its weight also decreases dramat- 
ically and it becomes difficult to see. The second panel 
of Fig. ETl shows A{k = tt, uj) vs. uj, on a linear scale. At 
A = 0, there is only one peak at +2t, as expected. As 
the coupling is turned on, this weight seems to spread 
around in a rather featureless, broad continuum. In fact, 
on a logarithmic scale (panel (c) of Fig. I^Tj) one can 
see that k = n weight is pulled down into all the bound 
states, in agreement with the previous data we showed. 
This weight, however, is so small that it is not visible 
on the linear scale. Finally, the lowest panel of Fig. [?T1 
shows the total spectral weight, or DOS. At A = we see 
the usual ID DOS, with the singularities near the band- 
edges rounded off because of the finite r] used. As A is 
turned on, one can recognize both the contribution from 
the k = and k = ir states to the total DOS: each bound 
state has a finite bandwidth due to its dispersion (this is 
to be contrasted to the upper panels, where the bound 
states are true delta-functions, with a width defined by 
the broadening 77). As the coupling strength increases, 
the number of bound states increases; they are spaced 
by roughly ft, their bandwidths narrow down, and their 
weights also decrease. In other words, they approach the 
expected Lang-Firsov behavior. 

Thus, this figure actually answers the question posed in 
the Introduction, regarding the evolution of the spectral 
weight from that of a free electron towards that of the 
impurity limit. While the MA results are certainly not 
exact, we can claim with a high degree of certainty that 
the main features are accurately captured, especially in 
the weak and in the strong coupling limits. This suffices 
to understand the physics of this problem, and given the 
simplicity of the approximation, to investigate in detail 
any number of other quantities we have not shown here, 
such as the self-energy. Of course, if one is interested 
in exact results for some particular set of parameters. 
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FIG. 22: (color online) 2D results for Q = 0.5t and ry = 0.04t. 
Top: A(k = (0, 0),(jj) vs. u and A. Bottom: total 2D spectral 
weight A((jtj) \s. u and A. The intensity scales are shown to 
the right of each plot. 



numerically-intensive methods have to be used. 

Qualitatively similar plots are obtained in 2D and 3D, 
as shown in Figs. [23 and ESI Of course, the A = DOS 
is very different, with a van-Hove singularity dit u = 
for the 2D case, and the characteristic nearest-neighbor 
hopping DOS in the 3D case. However, the appearance 
of multiple bands below the continuum as A increases and 
all the remaining phenomenology is very similar. Thus, 
we see no evidence of any qualitative differences in the 
polaron physics due to different dimensionality. 



V. SUMMARY AND CONCLUSIONS 

In this paper we analyzed the Green's function of the 
Holstein polaron, using the momentum average approx- 
imation, which consists in summing all the self-energy 
diagrams, but with each individually averaged over all 
its free propagator momenta. The resulting self-energy 
can be written as an infinite continuous fraction that is 
numerically trivial to evaluate. This procedure becomes 
exact in the limits ^ = and t = 0. 

We gauged the accuracy of this approximation by com- 
puting its corresponding spectral weight sum rules and 
comparing them against the exact sum rules, which are 
known for this type of Hamiltonian. We showed that the 
MA spectral weight satisfies exactly the first 6 sum rules. 
Even though this is quite impressive at first sight, it is 
actually no guarantee of overall accuracy, as the case of 
SCBA demonstrates. The SCBA spectral weight always 
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FIG. 23: (color online) 3D results for Q = 0.5t. Top: A(k = 
(0,0,0),cj) vs. 00 and A, for ry = 0.15t. Bottom: total 3D 
spectral weight A{co) vs. co and A, for 77 = 0.04t. The intensity 
scales are shown to the right of each plot. 

satisfies exactly the first 4 sum rules, even at large cou- 
plings A where it predicts very wrong results! The mean- 
ingful criterion of accuracy for the sum rules is to show 
that the vast majority of terms in all sum rules ^ and in 
particular the dominant terms in the various asymptotic 
limits, are captured by the approximation. MA indeed 
satisfies this very restrictive criterion. 

The accuracy of the approximation was also tested by 
direct comparison with data obtained through numeri- 
cally intensive methods. In all cases we obtain remark- 
able agreement, especially considering the ease of eval- 
uation of the MA results. The MA approximation is 
not exact and some features are not correctly captured 
(for example, the continuum starting at £^0 + ^) how- 
ever, in all cases, all the higher- weight features in the 
spectral weight are quantitatively and qualitatively well 
described by the MA approximation. Trading some of 
the accuracy of numerically exact but time consuming 
methods in exchange for very fast results which capture 
the main features accurately is a useful approach when 
trying to understand the main aspects of the physics of a 
problem, as well as when one is concerned about compar- 
ison with experiments. It is very unlikely that ARPES 
could capture very low- weight features, these would be 
lost in noise statistics. Thus, an approximation like MA, 
that quickly but accurately estimates results is very use- 
ful, to be followed, of course, by detailed numerics for the 
parameter sets of interest. 

It is important to note that the existence and accu- 
racy of approximations like MA is not guaranteed, in 
fact it can be regarded as a surprise for the case of the 
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Holstein polaron, that has been under investigation for 
almost five decades. However, this demonstration of its 
existence and efficiency in the Holstein polaron problem 
gives some hope for making progress for the general class 
of strongly-correlated systems problems. One can always 
use some ffavor of perturbation theory to understand be- 
havior in asymptotic limits, but the really challenging 
problems are set in regimes where perturbation does not 
apply. The MA approximation suggests that one way to 
make no n- trivial progress is to sum all diagrams, with 
each simplified enough so as to make the calculation fea- 
sible, but not so much as to really alter the physics. This 
is a very different approach from the usually employed 
summation of a subclass of diagrams. Note that there 
are many classes of problems with diagrams similar to 
the ones arising in the single polaron problem, although 
of course with different propagators and/or vertices. 

A first possible generalization of this work is to mod- 
els with several phonon branches, and/or momentum- 
dependent coupling ^q, and/or dispersive phonons, l^q. 
The way to achieve this for a single polaron is briefiy 
discussed in Ref. [s^- Given the length of this article, 
we postpone this discussion for future publications where 
results can also be shown. Other directions of generaliza- 
tion are to finite particle densities and/or finite temper- 
atures, and indeed to Hamiltonians which also include 
electron-electron interactions. 

Such work is currently in progress. It is still far from 
clear which cases will admit useful generalizations, how- 
ever the proof of existence of this method for the Holstein 
polaron problem is in itself encouraging. 
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APPENDIX A: COMPARISON WITH DMFT 

The MA self-energy of Eq. (|T^ looks similar to the 
DMFT self-energy, discussed in Ref. 32. This is not so 
surprising, since both become equal to the exact Lang- 
Firsov limit if the bandwidth goes to zero, and this can 
be re-written as an infinite continued fraction, as shown 
in Eq. ([T^ . However, this similarity may raise questions 
about the relationship between the two approximations. 
In this appendix, we briefly show that the two approxi- 
mations are very different at all finite t and g. 

The meaning of the MA go{^ and of corresponding 
DMFT Go{uj) (notation of Ref. Isl) is very different. The 
DMFT Go{uj) is obtained by solving exactly the problem 
of an impurity coupled to an environment in the d ^ oo 
limit, and then imposing the self-consistency condition 
that the impurity site behaves similarly to all other sites 
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FIG. 24: (color online) Comparison between the function 
Go (to) entering the DMFT self-energy, for cj/t = 0.5 and 
A = 0.5 (green) and A = 1.5 (yellow), and ^o(^) entering 
the MA self-energy, for the d ^ oo semi-elliptical DOS. 



in the environment. The DMFT Go(^) is calculated self- 
consistent ly using the following steps »^ (i) with some 
initial guess for Go{uj)^ one calculates the DMFT self- 
energy, given by a formula similar to Eq. OHJ^ with 
^o(^) replaced by Go{u); this self-energy is then used 
to calculate the total Green's function 
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where po(e) is the density of states of the non-interacting 
electrons. The usual procedure is to take d ^ oo and 
use as DOS the semielliptical formula corresponding to 
an infinitely branched Bet he tree. 



where W is the bandwidth for the non-interacting sys- 
tem. Then, (iii) the new Go(^) is extracted from the con- 
dition that G{uj) = [Gq^{uj) — S(a;)] and the proce- 
dure is repeated until self-consistency is reached. Reach- 
ing self-consistency is a non-trivial numerical task, espe- 
cially compared to obtaining the MA ^o(^) (see below). 
More importantly, the DMFT Go(^) is an explicit func- 
tion of g and Q. 

In contrast, the MA ^o(^) is the momentum average 
of the free propagator: thus, it is a known function, in- 
dependent of the parameters g and In particular, for 
the semi-elliptical DOS, we have simply: 
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A comparison of these functions is provided in Fig. 
They are clearly different. Moreover, note that in 
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the MA approximation, G(k, uj) is an explicit function of 
the momentum. The MA self-energy for the Holstein po- 
laron happens to be independent of the momentum, but 
this a consequence of the simplicity of the model, not an 



"in-built" feature like in DMFT. Generalizations to mod- 
els which have a momentum-dependent coupling and/or 
dispersive phonons, lead to momentum-dependent self- 
energies^ 
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